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Abstract 

This paper describes an AERONET-based Surface Reflectance Validation Network (ASRVN) 
and its dataset of spectral surface bidirectional reflectance and albedo based on MODIS TERRA 
and AQUA data. The ASRVN is an operational data collection and processing system. It 
receives 50x50 km 2 subsets of MODIS LIB data from MOD APS and AERONET aerosol and 
water vapor information. Then it perfonns an accurate atmospheric correction for about 100 
AERONET sites based on accurate radiative transfer theory with high quality control of the input 
data. The ASRVN processing software consists of LIB data gridding algorithm, a new cloud 
mask algorithm based on a time series analysis, and an atmospheric correction algorithm. The 
atmospheric correction is achieved by fitting the MODIS top of atmosphere measurements, 
accumulated for 16-day interval, with theoretical reflectance parameterized in terms of 
coefficients of the LSRT BRF model. The ASRVN takes several steps to ensure high quality of 
results: 1) cloud mask algorithm filters opaque clouds; 2) an aerosol filter has been developed to 
filter residual semi-transparent and sub-pixel clouds, as well as cases with high inhomogeneity of 
aerosols in the processing area; 3) imposing requirement of consistency of the new solution with 
previously retrieved BRF and albedo; 4) rapid adjustment of the 16-day retrieval to the surface 
changes using the last day of measurements; and 5) development of seasonal back-up spectral 
BRF database to increase data coverage. The ASRVN provides a gapless or near-gapless 
coverage for the processing area. The gaps, caused by clouds, are filled most naturally with the 
latest solution for a given pixels. The ASRVN products include three parameters of LSRT model 
(, k L , k G , k v ), surface albedo, NBRF (a normalized BRF computed for a standard viewing 
geometry, VZA=0°, SZA=45°), and IBRF (instantaneous, or one angle, BRF value derived from 
the last day of MODIS measurement for specific viewing geometry) for MODIS 500m bands 1- 
7. The results are produced daily at resolution of 1 km in gridded format. We also provide cloud 
mask, quality flag and a browse bitmap image. The new dataset can be used for a wide range of 
applications including validation analysis and science research. 

1. Introduction 

Validation of the moderate resolution (~lkm) surface reflectance products, including 
spectral bi-directional reflectance factor (BRF), albedo, vegetation indices and others, is an 
important component of the Earth Observing System (EOS) [1] and the National Polar Orbiting 
Environmental Satellite System programs. Its goal is to establish the accuracy of environmental 
data products on regional and global scales for a broad range of atmospheric and surface 
conditions. The EOS program has developed a multi-level strategy with a strong field campaign 
component [2-5]. The field measurements required for direct validation analysis provide a 
detailed and comprehensive look at the local properties, but they usually involve significant 
resources, are subject to weather uncertainties, and are strongly limited in temporal and spatial 
coverage. Due to these constraints, recent validation efforts have proposed that product accuracy 



assessment should also utilized a globally representative sample of sites to complement the direct 
validation sites [6], This concept has been endorsed by the Committee on Earth Observing 
Satellites (CEOS) as the Bench-mark Land Multisite Analysis and Intercomparision of Products 
(BELMANIP) [7], 

In this paper, we present an alternative validation approach for moderate resolution 
global surface reflectance products over Aerosol Robotic Network (AERONET) sunphotometer 
sites [8], The idea is to collect the best ancillary infonnation on atmospheric aerosol and water 
vapor, and perform an independent atmospheric correction of satellite measurements based on 
accurate radiative transfer theory with high quality control of the input data and results. In the 
past several years, we have implemented this idea in the AERONET-based Surface Reflectance 
Validation Network (ASRVN) which is an automated data collection and processing system 
residing on a dedicated workstation. ASRVN operationally receives the satellite sensors’ LIB 
data (currently MODIS TERRA and AQUA from the Goddard’s MOD APS and MISR from 
Langley’s DAAC) and aerosol and water vapor information from the AERONET server. After a 
successful test of data integrity and completeness, ASRVN automatically performs rigorous 
atmospheric correction on each sensor’s data, creating a sensor-specific record of spectral BRF, 
albedo and derivative products (e.g., vegetation index) over more than 100 AERONET sites 
globally. For each site, products are stored in a gridded format at 1 km resolution for an area of 
50x50 km 2 for MODIS and 32x32 km 2 for MISR. Currently, the ASRVN dataset contains 8 
years of MODIS TERRA data (since 2000), and 4 years of MODIS AQUA data. The AQUA 
dataset will be completed together with MODIS Collection 5 land re-processing, which produces 
the MODIS AQUA subsets for the ASRVN. 

Many of the ideas implemented in ASRVN originated in the MODIS BRDF/Albedo 
algorithm [9] similarly based on the time series analysis within 16-day intervals. In this paper, 
we will be using the term BRF (bi-directional reflectance factor) rather than more commonly 
used BRDF (bi-directional reflectance distribution function), which is just a factor of n smaller 
than BRF. 

The ASRVN algorithm for MISR data was presented previously [10]. This paper 
describes processing algorithm for MODIS measurements (sections 2-3). Examples of ASRVN 
dataset are given in section 4. The paper concludes with a summary. 


2. ASRVN Infrastructure 

The BRF retrievals from MODIS data use several clear-skies measurements acquired on 
successive observations of the same area under different viewing geometries. The ASRVN 
processing starts by gridding the latest MODIS swath data (generated operationally over each 
AERONET site by the MOD APS production system) and placing the resulting ASRVN Tile in a 
16-day processing Queue, which implements a sliding temporal window algorithm. Note the 
ASRVN Tile is distinct from the much larger MODIS Land Team’s production Tile. Next, it 
finds a relevant AERONET aerosol and water vapor record and computes radiative transfer 
functions required for atmospheric correction. After ensuring the quality of input data by filtering 
clouds and spatially variable aerosols, the ASRVN algorithm performs atmospheric correction 
and checks the quality of the final solution. The next three sections will describe pre-processing 
steps preceding the BRF retrieval. 



2.1 Ancillary Data 

AERONET sunphotometers sample the direct solar radiation each 15 minutes, and 
sample diffuse sky radiance over a wide range of angles every 60 minutes during the daytime. 
AERONET ’s automated processing system generates AOT and column water vapor from the 
direct solar measurements. Typical AOT uncertainty for a field instrument is 0.01-0.02 and is 
spectrally dependent. The inversion algorithm [11] uses almucantar sky measurements to retrieve 
aerosol microphysical properties (particle size distribution and refractive index) and 
concentration. AERONET applies several tests to ensure reliability of retrievals, such as SZA > 
45°, AOTo. 44 > 0.4, and that there were at least 21 independent angles used in each inversion. The 
tests analyze sensitivity of retrievals to the single scattering albedo, and to the phase function at 
large scattering angles. These quality assurance tests significantly reduce the number of complete 
aerosol characterization records, as compared to the number of AOT records. 

The ASRVN algorithm starts with the selection of AERONET aerosol optical thickness 
and column water vapor values within 30 minutes of satellite overpass. If these conditions are 
met, the algorithm selects the inversion record with aerosol microphysical parameters and size 
distribution within 2 hours of overpass. Otherwise, it uses an aerosol climatology model for a 
given location derived from multi-year AERONET statistics of reliable retrievals [12]. Because 
full AERONET inversions are less accurate when aerosol concentration is low, the climatology 
background aerosol model is always used in our algorithm for clear atmospheric conditions 
(currently defined as AOTo,44<0.3). Our testing demonstrates that the aerosol climatology 
significantly improves the stability of the time series of derived surface albedo. 

Following the selection of aerosol parameters, ASRVN algorithm calculates the aerosol 
optical thickness (AOT M ie), single scattering albedo, and scattering phase function using a look- 
up table approach [13]. Depending on AERONET sphericity index, either a spherical (Mie) 
aerosol model or a model of spheroids is used. The above calculations provide the spectral 
dependence of extinction (AOTmie) in the MODIS wavelengths. However, the AOT from direct 
solar measurements may differ from AOTmie for a number of reasons, from the time difference 
between inversion and direct AOT measurement to uncertainties associated with the AERONET 
inversion algorithm [14] etc. For this reason, AOTmie is further scaled by fitting it to the 
measured AOT at three AERONET wavelengths (0.44, 0.67, and 0.87 pm). Once aerosol optical 
parameters are defined, the radiative transfer model SHARM [15] calculates the required 
radiative transfer functions for the specific water vapor and spectral response functions of 
MODIS TERRA or AQUA instrument using the Interpolation and Profile Correction (IPC) 
method [16]. 

2.2 Implementation of Time Series Processing 

To execute time series processing (sliding window algorithm), ASRVN first grids MODIS 
level IB (LIB) calibrated and geolocated data to a regular 1 km grid. We use the MODIS land 
gridding algorithm [17] with minor modifications that allow us to better preserve the angular 
anisotropy of signals in the gridded data when measured reflectance is high, for example over 
snow, thick clouds or water with glint. Next, gridded MODIS data (Tiles) are placed in the 
processing Queue, which can hold up to 16 days of successive measurements. The ASRVN 
processing uses both individual grid cells, also called pixels below, and fixed-size (25x25 km ) 
areas, or blocks, required by the cloud mask algorithm. In order to organize such processing, we 
developed a framework of C++ classes and structures (algorithm-specific Containers). The class 



functions are designed to handle processing in the various time-space scales, for example at the 
pixel- or block-level, and for a single (last) day of measurements or all available days in the 
Queue, or for a subset of days which satisfy certain requirements (filters). The data storage in the 
Queue is efficiently organized using pointers, which avoids physically moving the previous data 
in memory when the new data arrive. 

The structure of the Queue is shown schematically in Figure 1. For every day of 
observations, MODIS measurements are stored as Layers for reflective bands 1-7 and thermal 
band 32, all of which are required by the CM algorithm. Besides storing gridded MODIS data 
{Tiles), the Queue has a dedicated memory (Q-memory) which accumulates ancillary 
information about every block and pixel of the surface for the cloud mask algorithm {Refcm data 
structure). It also keeps information related to the history of previous retrievals, such as surface 
BRF parameters and albedo. Given the daily rate of MODIS observations, the land surface is 
relatively static over most 16-day periods. Therefore, knowledge of the previous surface state 
significantly enhances both the accuracy of the cloud detection and the quality of atmospheric 
correction by imposing a requirement of consistency of the time series of BRF and albedo. 

2.3 Data Quality Control: Cloud Mask and Aerosol Filter 

From the start, ASRVN was designed to work with a multipixel area rather than a single 
pixel centered at an AERONET location in order to provide optimal visual control over the input 
and output results. Visual analysis of RGB images is superior for complex data quality 
assessment and troubleshooting situations which is rarely achievable with a pixel- (or point-) 
level analysis. 

Although AERONET produces an internal cloud mask [18], it is not sufficient for 
atmospheric correction over the 50x50 km processing area. For example, the sunphotometer 
may provide AOT measurements from a direct sun view through a gap in the clouds. Usually, the 
sunphotometer’ s time of measurement differs from satellite overpass time, in which case changes 
in cloudiness may have occurred. For these reasons, we implemented a new MAIAC Cloud Mask 
algorithm [19] as part of the ASRVN processing. This algorithm uses the time series analysis and 
an image-based rather than pixel-based analysis. It explicitly addresses the problem of cloud 
searching by identifying a clear-skies comparison target. MAIAC CM constructs the reference 
clear-skies image of the surface and stores it along with other ancillary information about 
reflectance and brightness temperature for every surface block (25x25 km ). This ancillary 
information, required for cloud detection, is continuously updated with the latest cloud-free 
measurements. This allows CM to dynamically adjust to gradual (seasonal) and rapid surface 
changes caused by snowfall, fires, floods, etc. The CM algorithm is enhanced by an internal 
dynamic land-water-snow classification, which allows processing flexibility over varying surface 
types. 

Our experience with MODIS data processing dictated the need for additional data 
screening. This is required when aerosols have significant spatial variation and a single 
AERONET AOT value does not represent the full processing area, as well as in cases of 
undetected, usually semi-transparent and sub-pixel, clouds. This screening was implemented 
through an “aerosol filter”. Using the known surface BRF from previous retrievals, the algorithm 
computes the pixel-level AOT in the blue band from the latest MODIS measurements. The AOT 
retrieval is a fast algorithm based on a look-up table pre-calculated for a standard continental 
aerosol model. The computed AOT is used solely to assess spatial homogeneity of aerosols over 



the processing area, and to find deviations, which usually indicate previously undetected clouds 
and sometimes spatially varying aerosols. Specifically, the algorithm generates an AOT 
histogram from the non-cloudy pixels, filters the highest 20% values as possibly cloudy, and 
finds the average value (AOT av ) for the remaining 80% of the pixels. This AOT av is assumed to 
represent the average clear-skies aerosol loading over the processing area, which should 
correspond to the AERONET AOT. The AOT av is used next to further filter “suspicious” 
(contaminated) data as follows: if the atmosphere is clear (AOT av <0.25), then the algorithm 
filters only pixels with the high AOT values exceeding the average by 0.15 or more. This 
threshold was determined through trial and error. Otherwise, it filters the high and low values 
symmetrically if the difference with AOT av exceeds ±0.15. This relatively simple technique 
allows us to filter sub-pixel clouds, contrails and other fonns of thin cirrus and semi-transparent 
clouds. With the introduction of this additional filter, we witnessed a dramatic enhancement in 
the quality of th e ASRVN atmospheric correction. 

Figure 2 shows examples of the cloud mask algorithm and aerosol filter for the Goddard 
Space Flight Center (Greenbelt, Maryland, USA) site. It shows that the CM algorithm captures 
most of the opaque clouds, whereas the aerosol filter captured additional sub-pixel and semi- 
transparent clouds and cases of spatially variable aerosols. The example at the bottom of Figure 2 
also shows the detection of cloud shadows by MAIAC CM algorithm. Shadows are detected with 
a simple threshold algorithm which compares the latest MODIS measurement (p meas ) with 
predicted reflectance from the previously retrieved BRF model parameters (p preel ): 

IF p meas < p pred - 0.12 => CLOUD SHADOW. 

Here, we use MODIS wavelength 1.24 pm (band 5) which experiences minimal atmospheric 
distortions and is usually bright over land so the change of reflectance due to cloud shadow can 
be easily detected well above the sensor noise level. 

3. Atmospheric Correction 

Once the cloud mask, enhanced by aerosol filter, is applied, the ASRVN algorithm filters 
the time series of MODIS measurements for every single pixel, and places the remaining good 
data in a “container”. The container stores measurements along with computed RT functions for 
the cloud-free days of the Queue. If the number of good measurements exceeds 3 for a given 
pixel (see section 3.2), then the coefficients of Li-Sparse Ross-Thick (LSRT) BRF model [20] 
are computed. The LSRT model is used in the MODIS BRDF/albedo algorithm [9]. 

3.1 Inversion for LSRT Coefficients 

In the operational MODIS land processing system, the BRF is detennined in two steps: 
first, the atmospheric correction algorithm derives surface reflectance for a given observation 
geometry using a Lambertian approximation [21]. Next, three LSRT coefficients are retrieved 
from the time series of surface reflectance accumulated for a 16-day period [9]. The Lambertian 
assumption simplifies the atmospheric correction but imparts biases which depend on 
observation geometry and atmospheric opacity. Tests show that the Lambertian assumption leads 
to a more Lambertian BRF shape while the true BRF is more anisotropic [22]. 

The ASRVN algorithm derives surface LSRT coefficients directly by fitting the radiative 
transfer solution to the measured TOA reflectance accumulated over a 4-16 day period. The 
inversion is based on a high accuracy semi-analytical Green’s function solution [23, 24], which 



in combination with LSRT BRF model provides an explicit parameterization of TO A reflectance 
in terms of the surface BRF model parameters K = {k L , k°, k 1 j 7 . According to the derivation 
provided in the Appendix, the TOA reflectance can be expressed as: 

R(ju 0 , ju,(p) = R D (ju 0 , ju,(p) + k L F L (ju 0 , ju) + k G F G (ju 0 , ju,(p) + k 1 F* (ju 0 , ju,(p) + R nl (ju 0 , ju) , ( 1 ) 

where R° is the atmospheric (path) reflectance, and R" 1 is a small non-linear term proportional 
to the product of surface albedo and spherical albedo of the atmosphere ( R" 1 oc qc 0 ). Functions 
F l , F 1 and F G depend on geometry and atmospheric conditions. They are weakly non-linear 
in ^'-coefficients through the multiple reflection factor a = (1 - qc o y\ 

The quasi-linear form of equation (1) leads to a very efficient iterative minimization algorithm: 

RMSE= Y (tj n) - F G k L(n) - F; k V(n) - F G k G(n) ) 2 =min, r {n) =R-R° (2) 

J J J jJ {K} 

where index j denotes measurements for different days, and n is the iteration number. Equation 
(2) provides an explicit least-squares solution for the kernel weights. In matrix form, the solution 
is written as: 

K (n) = A-'b {n) , (3) 

where 
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In the first iteration, the small non-linear term is set to zero, R" H0) = 0 , and the multiple 

reflection factor a (see Appendix) is set to one, a f0) = 1 . These parameters are updated once 
after the BRF coefficients are calculated in the first iteration. Except for snow-covered surfaces, 
the problem converges with high accuracy in two iterations because the non-linear terms are 
small. Currently, the ASRVN algorithm does not make retrievals over snow. 

The described algorithm has a high computational efficiency. Compared to the radiative 
transfer computations, the time required to evaluate functions F m ( m=L , V, G) and R" 1 is 
negligible. The integrals required for these functions (see Appendix) need to be calculated only 
once regardless of the number of iterations. Finally, small variations of the viewing geometry 
across the processing area are neglected and the RT calculation is done once per observation 
using the geometry of the central pixel. 


3.2 Solution Selection and Update 

Although the LSRT model leads to an efficient BRF retrieval algorithm, there are several 
caveats associated with this model. The LSRT kernels are not orthogonal, are not positive-only 
functions, and are normalized in a somewhat arbitrary fashion that is not linked to radiative 
transfer theory. These factors reduce the stability and uniqueness of the solutions, such that small 



perturbations in measurements may lead to significantly different solutions. The high goodness- 
of-fit at the measurement angles does not guarantee the correct shape of the retrieved BRF, and 
may result in negative BRF values at other angles. The albedo, being an integral function of 
BRF, is especially sensitive to an incorrect BRF shape. For these reasons, we developed several 
tests to remove unrealistic solutions. 

The initial validation of the solution (see Figure 3) checks that the maximal difference 
over all days of the Queue between measured and computed TO A reflectance does not exceed a 
specified threshold ( | R Meas - /? /v ' y | >0.08). The day (measurement) with the highest deviation is 

excluded from the Queue and the inversion is repeated. If the number of measurements goes 
below three after the exclusion, no retrieval will be made for this pixel. 

If a solution provides a good agreement with measurements for all days, the algorithm 
verifies that values of the direct-beam albedo ( q ) at SZA=15°, 45°, 60° are positive. Finally, the 

new solution must be consistent with the previous solution: |^(45°) - q Prev (45° )| <A( a), where A 

is the band-dependent threshold currently equal to 0.04 (blue), 0.05 (green and red), and 0.1 
(NIR and shortwave infrared bands). Consistency of the time series of BRF and albedo is 
characterized by a status index. Initially, the confidence in the solution is low {status= 0). Each 
time the new retrieval agrees with the previous retrieval, status increases by 1. When status> 3, 
the retrieval is considered reliable. 

When the new solution is validated, the coefficients of the BRF model and direct-beam 
albedo <y(45°), stored in the Q-memory, are updated. The update is done with relaxation, 
designed to mitigate random noise of retrievals: 

kf w = ( kf w + k Prev )/2. (4) 

This method of update increases the quality of the BRF and albedo product when the surface is 
relatively stable, but it delays the response of the solution to surface changes. 

Often, the solution for some pixels or the full area cannot be produced because of lack of 
clear-skies measurements. In these cases, we assume that the surface does not change and we 
fill-in data gaps with the previous solution for up to a 32-day period. In most cases, this 
assumption of a stable surface is reasonable. The gap-filled pixel is marked as “Extended” in the 
quality assurance (QA) value, with parameter QA .nDelay giving the number of days since the 
last update (see section 3.6). 

3.3 ASR VN Products 

The ASRVN computes two main products at 1 km resolution for seven 500m MODIS 
bands, the set of BRF coefficients, and the surface albedo. The albedo is defined by Equation (A- 
5a) as a ratio of surface-reflected to incident radiative fluxes. Thus, it represents a true albedo at 
a given solar zenith angle in ambient atmospheric conditions, the value, which can be directly 
compared to ground-based measurements. 

ASRVN also computes several derivative products useful for science data analysis and 
validation: 

1) NBRF - a BRF Normalized to the common geometry of nadir view and SZA = 45°. 
This product is analogous to MODIS NBAR (nadir BRF-adjusted reflectance) product (part of 



the MOD43 standard product suite). With the geometry variations removed, the time series of 
NBRF is useful for studying vegetation phenology, perfonning surface classification, etc. 

2) IBRF - an instantaneous (or one-angle) BRF value for the specific viewing geometry 
of the last day of observations. In essence, IBRF is a reflectance which would be measured if the 
atmosphere were absent. This product is calculated from the latest MOD1S measurement 
assuming that the shape of BRF, known from previous retrievals, didn’t change. To illustrate 
computation of IBRF, we re-write equation (1) for the measured TO A reflectance as follows: 

R(jU 0 ,n,(p) = R D (ii 0 ,jU,(p) + bR Surf ( fi () , fi, (p) , (5) 

where R Swl combines all surface related terms and can be calculated using the previous solution 
for BRF (BRFX) and AERONET aerosol infonnation. b is spectrally-dependent scaling factor. 
Then, 

IBRF \ (/r 0 ,n,<P) = b.BRF, (// 0 , ju, <p) . (6) 


This algorithm (Equation 6) will be referred to as scaling. This description was given for the 
purpose of illustration. In reality, R Sl " j is a non-linear function so computing parameter b)_ and 
IBRF is done accurately using the formulas given in the Appendix. 

The algorithm computing scaling coefficient (and IBRF) is shown in Figure 3 on the 
right. First, the algorithm filters measurements which differ from the theoretically predicted TO A 
reflectance based on the previous solution ( R' ( f RT ) by more than factor of A(X). Then, the scaling 

coefficients are computed, and the consistency requirement is verified as: 0.8<Z? / i<l .2. If all 
conditions are satisfied and the status of the pixel is high (status> 3), then the LSRT BRF 
parameters of the pixel stored in the Q-memory are updated with the scaled solution: 
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Based on its definition, IBRF is well suited for validation of surface reflectance product of the 
operational atmospheric correction algorithm (standard product MOD09 [21]). 

The cloud mask, the RGB browse images for MODIS TOA reflectance and for the 
ASRVN NBRF and IBRF, and the QA flag, are also standard parts of the ASRVN product suite. 
The products are saved in HDF-EOS format files, which can automatically keep geolocation 
information and allow data to be ported to virtually any computer platfonn, regardless of the byte 
order of the native platform. A list of the ASRVN products and their MODIS products 
counterparts are given in Table 1, and the definition of the cloud mask values is provided in 
Table 2. 


3.4 Update in Case of Rapid Surface Change 

Time series processing is intrinsically controversial when a surface changes rapidly. 
Since inversions are generally ill-posed problems, one desires all available cloud-free 
measurements and a maximal time window in order to reduce the RMSE. This approach, which 
reduces the impact of noise in the data (including that of gridding and of residual clouds) and 
which ensures a more robust BRF shape, is best when the surface is stable throughout the 
accumulation period. For example, this is the case for natural ecosystems in mid-latitudinal 
summers. In contrast, detecting and tracking surface changes like spring green-up, agricultural 



harvesting or fall senescence requires the least possible number of days in the inversion Queue. 
Such retrievals may have a considerable amount of spatial and spectral noise. Indeed, it is 
difficult to assess the reliability of solutions when the surface reflectance is changing rapidly or 
abruptly, especially given the possibility of data gaps due to clouds. 

While the response of the 16-day solution may be delayed, the IBRF tracks spectral 
changes immediately. The update of Q-memory with the latest measurements by equation (7) 
was found to significantly accelerate the response of LSRT coefficients ( K - ), and hence of the 

NBRF, to changing surface conditions. Yet, in some cases we found this insufficient. To amend 
results in these cases, we added a separate update of the Q-solution with IBRF based on change 
detection. This path is shown at the bottom of Figure 3. Usually, seasonal surface changes related 
to the vegetation cycle are accompanied by a correlative change in the red-NIR bands, for 
example, a simultaneous decrease of red and increase of NIR reflectance during green-up, and 
the opposite changes during fall senescence and defoliation for broadleaf forests. In our case, the 
change detection is based on the top-of-canopy Normalized Difference Vegetation Index 
(NDVI). The scene-average NDVI, normalized to a standard viewing geometry of NBRF, is 
calculated using the IBRF (NDVI 1 ). It is compared to the scene-average NDVI Q calculated using 
the NBRF stored in the Q-memory (previous reliable solution). A change is defined as when the 
difference between the two values exceeds ±0.01. At the same time, the red and NIR reflectance 
should change accordingly. For example, the following set of rules defines the change during 
spring green-up: 

NDVI 1 - NDVI 0 > 0.01 , - PL < 0 . p' m - Pj?„ > 0 . 

Such an approach filters fluctuations in NDVI caused by factors unrelated to the vegetation 
signal, such as residual cloudiness or undetected cloud shadows. Once the change is detected, the 
Q-memory is updated using Eq. (7). Empirically, we found that this approach, which is 
responsible for 10-15% of all updates, performs robustly at different global AERONET 
locations, for example in North America, Eurasia, or Africa. 

The NBRF and albedo are also updated following an update of K x coefficients. This method, 

which requires only one last day of measurements when the BRF shape is reliably known, is 
similar to the scheme used in the MODIS BRDF/albedo algorithm [9] and its reliability has been 
tested with a global composite AVHRR data set [25]. An advantage of this strategy is that once 
initialized, the algorithm provides near gapless coverage, and a continuous time series of NBRF 
and albedo. 

Besides providing a faster response to surface change, our update strategy assures fast 
removal of retrieval artifacts, mainly residual clouds, which is the most common problem. 

3.5 Seasonal Ancillary BRF 

When no reliable retrievals are made during the past 32 days, which is usually caused by 
high cloudiness, the previous retrievals are considered unreliable and the Q-memory is refreshed 
with fill values. After that, it may take the algorithm a considerable amount of time to re- 
initialize, during which time no results will be produced. To remedy this situation, we developed 
the historical 1 km resolution BRF database from 5 year retrievals, which is used for BRF- 
scaling when the Q-memory is being refreshed. The database contains one set of spectral K- 
coefficients for every pixel for each of the four calendar seasons. Initially, the database is built 



from a multi-year run of ASRVN. The images used for averaging, are selected according to the 
standard deviation of NBRF (onbrf) computed for the processing area. Empirically, we found 
that with our data quality assurance, the top 20% of images with the highest standard deviation 
usually contain artifacts from residual clouds or unreliable BRF solutions, whereas the remaining 
80% of retrievals have a good quality. Thus, for each site and every season, we first generate the 
histogram of ctnbrf, find the 80% threshold, and average the images with a-values lower than the 
threshold. 

Once the seasonal BRF database is created, it is supported by an offline background 
algorithm which updates it with the latest good quality solution. 

3.6 Quality Assurance (QA) Flag 

For each execution, the algorithm creates a pixel-level QA flag to indicate the overall 
quality and the internal processing path. The QA information consists of 16-bit compound bit- 
fields, as summarized in Table 3. 

The QA. overall field indicates the overall data quality. Four values are possible: 

1) Good. This means that the status of solution is 3 or greater. In other words, at least three 
of the last LSRT retrievals agreed, and the calculated IBRF was found to be consistent 
with the model predictions; 

2) Acceptable. In this level, the status is not high (<=2), but the calculated IBRF agrees with 
the model prediction. This case may represent a good solution where only a few retrievals 
are available because of a gap in AERONET data or clouds; 

3) Extended. In this case, either the solution was not produced because of a lack of clear-sky 
measurements, or the calculated IBRF did not agree well with the model prediction. In 
this case, the previous reliable solution for a given pixel is used to fill-in the values of the 
LSRT coefficients and NBRF. The IBRF in this case is not produced; 

4) Not created. This usually happens at the beginning of processing when the Q-memory is 
not yet initialized. 

The field QA.nDelay gives the number of days since the last update of the Q-memory. If 
QA.nDelay =0, then this pixel contains the most recent retrieval. If the Q-memory was not 
updated for 32 or more days, the infonnation for a given pixel will be overwritten with the fill 
value. If the solution is calculated with (Eq. 7), the QA.scale field will be set to 1 to indicate that 
the solution is “scaled” from the previous reliable retrieval. The value of the cloud mask for a 
given pixel is stored in the field QA.cloud. We also mask pixels which are adjacent to cloudy 
pixels, where greater errors in the atmospheric correction are expected. For such pixels, the 
QA.adjCloud field is set to 1 (the default value is 0). 


4. Processing Examples 

We selected three AERONET sites, Goddard Space Flight Center (GSFC, 38.9925°N, 
76.84° W), Mongu, Zambia (15.25°S, 23.15°E) and Solar Village, Saudi Arabia (24.9 1°N, 
46.4 1°E) to illustrate the ASRVN dataset. These sites have very different land cover types, 
atmospheric conditions and seasonal variations. 



• The GSFC site is located in a northern suburb of Washington, DC. It is a mixture of 
urban residential area, small deciduous broad leaf forest stands and small patches of 
agriculture cropland. 

• The Mongu site is located on the eastern side of Zambezi River in western Zambia. 
The western part of the area is a floodplain mostly covered with grasses, and the 
eastern part is mainly a sandy soil with sparse vegetation. Significant biomass 
burning happens in August and September. In the wet season (November-March), this 
area has frequent cloud cover. 

• The Solar Village site is a desert area, with very stable surface conditions. The 
dominant aerosol source is dust. 

Below, we provide several examples to demonstrate ASRVN products, data quality and potential 
applications. For this study, we used MOD APS Collection 5 data exclusively. 

4,1 NBRF time series 

The NBRF product has been corrected for both atmospheric effects and variations of 
view geometry. Thus, NBRF variations should be closely related to changes in surface 
conditions. Figure 4 shows the seasonal dynamics of the NBRF images for the three sites. The 
“true color” images are composed with equal weights from the red, green and blue spectral 
bands. The columns show the gridded TOA reflectance, NBRF and the new cloud mask (the CM 
legend is given in the figure’s caption). Seasonal changes are easy to observe in the first two 
columns. The vegetative cover reaches a maximum in July-August for the GSFC site, whereas at 
Mongu, green vegetation is most abundant in January-February, at the height of the wet season. 
Surface reflectance at the Solar Village site exhibits little variation throughout the year as is 
expected for a desert surface with little vegetation. 

To demonstrate the algorithm’s performance with different surface types, we selected two 
pixels for GSFC that differ in the amount of vegetation - a “green” pixel [pixel (16, 36)] in the 
middle of a small deciduous forest stand to the north-east of the center, and a relatively bright 
“urban” pixel [pixel (46, 3)] representing a typical residential area with mixed vegetation at the 
lower-left corner of the image.. For comparison, we also selected a “bright” pixel [pixel (10, 20)] 
in the desert region of Solar Village. The locations of these pixels are indicated by color circles 
in Fig. 4. The NBRF time series for these pixels are shown in Fig. 5. For the “green” pixel, the 
NBRF in the NIR band increases rapidly during springtime, while the red, green, and blue 
reflectances decrease. There is an interesting dynamic between the red and green signals. In the 
early spring, reflectance in the red channel is greater than in the green channel. This is typical of 
most soils. With the spring green-up, the red band reflectance decreases significantly due to 
chlorophyll absorption while the change in the green band reflectance is much smaller. During 
the autumn season, the bands change in the reverse direction as expected for senescing 
vegetation (Fig. 5a). The “urban” pixel shows a similar temporal trend, but with much smaller 
amplitude (Fig. 5b). The time series of the “bright” pixel at Solar Village does not show much 
variation throughout the year. The NBRF of band 7 and the NIR, red, green and blue bands 
remains around 0.5, 0.42, 0.35, 0.24, and 0.14 respectively . The variation is about +/- 0.03-0.05 
in each band. 

The data gap in Fig. 5 in year 2004 is due to the incompleteness of the MODIS dataset 
we acquired to date. The remaining gaps will be filled in upon completion of the MODIS land 



Collection 5 reprocessing, which also generates the MODIS subsets for the AERONET sites for 
the ASRVN. 


4.2 IBRF vs. NBRF 

As discussed in section 3.4, the NBRF, which is retrieved from 16 days of measurements, 
may have a delayed response to surface changes. In contrast„the IBRF, derived from the last day 
of measurements, tracks surface spectral changes immediately. This improvement in temporal 
sensitivity is sometimes achieved with compromised data quality. 

Figure 6 shows a comparison of NBRF and IBRF for the GSFC site for the spring green- 
up period (days 95-122) of 2005. There is a 7-day gap in the AERONET records due to 
cloudiness after day 101. As a consequence, the NBRF and IBRF images show a noticeable color 
difference on day 108, representing a delayed response of the NBRF to surface change. With the 
accumulation of six additional measurements, the color of the NBRF and IBRF images become 
consistent again on day 117. The lag in NBRF depends on the number of available clear-skies 
measurements and the rate of surface change. 

Note, however, that the brightness of the IBRF images in Fig. 6 vary with the view zenith 
angle. This occurs because the IBRF images are not normalized, whereas the NBRF images 
show a stable green-up trend as expected. This confirms that IBRF is more useful for detecting 
surface changes as well as effects of storms, fires, etc. The more stable NBRF is likely more 
suitable for process models, where the models may be sensitive to noise in input data fields. 

4.3 NDVI 

The Normalized Differential Vegetation Index (NDVI) is a commonly used parameter to 
characterize vegetation canopies. NDVI can be directly derived from the ASRVN BRF and 
albedo products. NDVI can be generated from different reflectance parameters, and its sensitivity 
to viewing geometry and atmospheric conditions can change accordingly. Figure 7 shows the 
NDVI time series for the “green” and “urban” pixels of the GSFC site as calculated from NBRF, 
surface albedo, IBRF (representing different forms of top-of-canopy NDVI) and from TOA 
reflectance. One can see that, in agreement with expectations, the variation of NDVI derived 
from NBRF and surface albedo is much smaller than that derived from IBRF and TOA 
reflectance. The NDVIs derived from the atmospherically corrected products (NBRF, albedo and 
IBRF) are also greater than the TOA NDVI. This is an expected since atmospheric correction 
tends to reduce the red band reflectance more than the NIR signal. Comparing the NDVIs 
derived from IBRF and NBRF, we find that the NBRF NDVI generally responds to the seasonal 
changes in a timely manner, except for a delay in a some cases which are usually related to a gap 
in the AERONET data. 

4.4 MODIS Terra vs. Aqua 

ASRVN creates a sensor-specific time series record of the surface reflectance over 
AERONET sites. Fig. 8 compares the NBRF time series for 1.5 years of the “bright pixel” of the 
Solar Village site. The NBRFs are very close to each other which suggests a good relative 
calibration between the MODIS Terra and Aqua instruments. Generally, cross-calibration of 
sensors flying in different orbits is a very difficult task. With the accumulation of a longer time 
record and sufficient global statistics, the ASRVN dataset may become a valuable source for the 
cross-calibration analysis of different sensors. 



5. Discussion 

To date, the ASRVN system has been used mainly as a prototype and to verify more 
manually-derived MODIS and MISR results reported elsewhere. In the near future, the National 
Polar Orbiting Enviromnental Satellite System (NPOESS) era will begin (~ 2013-2026), 
preceded immediately by the NPOESS Preparatory Project (NPP) mission (launch -2010) [26] , 
2008 ]. NPP data should initially overlap that of EOS Aqua as both are polar afternoon orbiters 
(1330 equator crossing). Plans are already underway to incorporate ASRVN into EOS-to- 
NPP/NPOESS data comparisons, as well as NPP/NPOESS product validation and long-term 
stability analysis. 

Under NPP/NPOESS, the Visible Infrared Imaging Radiometer Suite (VIIRS) will 
replace MODIS and NOAA’s Advanced Very High Resolution Radiometer (AVHRR) as the 
nation’s wide-swath multispectral sensor. The 22-band VIIRS provides most of the spectral 
measurements and capabilities afforded by MODIS. VIIRS bands have a resolution of 375 m or 
750 m, and therefore are highly compatible with the existing ASRVN framework. 

VIIRS will provide source data for several relevant Environmental Data Records (EDRs; 
similar to MODIS L2 swath products), including AOT, instantaneous broadband albedo, and 
Enhanced Vegetation Index (EVI). The latter is a top-of-canopy 3-band deterministic parameter 
similar to NDVI. The NPOESS Program has rigorously defined accuracy and precision error 
specifications for each EDR. Surface spectral surface reflectance will also be produced and 
archived, although the Program considers this an Intermediate Product (IP). Among other things, 
the IP designation does not carry accuracy specifications. 

NPOESS EDRs are primarily intended for operational users (e.g., weather forecasting), 
and thus ground processing will occur much more rapidly than for EOS. Most EDRs will be 
available within 30 minutes of observation. Under such conditions, maintaining constant product 
quality will be more challenging than for EOS. Unfortunately, because NPOESS is an 
operational program, validation resources will be significantly more limited than during the EOS 
era. Nevertheless, EDR validation is highly important, both for EDR users and to verify that 
EDRs meet contractual error specifications. Further, the complete NPOESS program will likely 
include three successive afternoon-overpass platforms (nominally launching in 2010 with NPP, 
2013 and 2018). Relatively seamless EDR performance from platform to platform is critical for 
many users. 

The NPOESS Land Validation Team, together with the prime NPOESS contractor, are 
responsible for assessing EDR accuracy. The Validation Team has identified AOT and Surface 
Reflectance as the highest priority land-related products for post-launch validation. Further, the 
Team has defined an initial strategy that strongly emphasizes the use of operational field network 
datasets (e.g., AERONET + ASRVN) over field campaigns. The Team is currently coordinating 
activities with the prime contractor, and we envision that, similar to MODIS, the operational 
processing center will generate EDR spatial subsets of relevant EDRs for ASRVN in near real 
time. Several targeted post-launch field campaigns will be undertaken to verify ASRVN results 
and document performance. NPOESS validation activities and resources are expected to 
decrease over time as the program focuses on long-term performance monitoring. ASRVN is 
expected to be critical resource throughout. 



6. Concluding Remarks 

This paper presented a new operational data collection and processing system ASRVN, 
which was initially designed for validation of the surface reflectance products. ASRVN collects 
MODIS and MISR LIB data for 50x50 km 2 areas for about 160 AERONET sites, and performs 
independent rigorous atmospheric correction using AERONET aerosol and column water vapor 
data. The atmospheric correction is achieved by fitting 16-day (multi-angle) sets of MODIS TOA 
measurements with theoretical reflectance accurately parameterized in terms of coefficients of 
the LSRT BRF model. The algorithm has a thorough data quality analysis component, including 
a cloud mask, an aerosol filter, and the control of the time series consistency of surface BRF and 
albedo. 

The algorithm is optimized in terms of noise reduction and its ability to track both 
seasonal and rapid surface change. During stable surface conditions, such as periods of 
maximum greenness during summertime in the northern hemisphere, a 16-day BRF retrieval 
gives a solution characterized by low noise. When the surface changes rapidly (e.g., agricultural 
harvesting), the algorithm gives greater weight to the last MODIS measurement of the period 
thus quickly adjusts to the surface change. 

We have also developed a back-up mechanism to cover extended periods of high 
cloudiness, when the number of cloud-free MODIS measurements in 16-day Queue is 
insufficient for the BRF retrieval. In this case, we use BRF scaling. This approach requires 
knowledge of the BRF shape which usually comes from the previous solution stored in the Q- 
memory. However, if there were no retrievals in the past 32 days, the algorithm substitutes the 
Q-solution with the historic seasonal BRF, which represents a quarterly-average BRF over years 
of retrievals. This method allows us to rapidly initialize retrievals after long periods of 
cloudiness and provide continuous gap-filled imagery of high quality. 

The ASRVN suite of products includes three parameters of the LSRT model ( k L , k G , k v ), 
surface spectral albedo, NBRF (a BRF value computed for a standard viewing geometry, 
VZA=0°, SZA=45°), and IBRF (a “scaled” BRF value derived from the last day of MODIS 
measurements at specific viewing geometry). All parameters are produced daily for seven 500 m 
MODIS bands at gridded 1 km resolution. We do not store vegetation indices like NDVI which 
are easy to produce from the available data. 

The ASRVN dataset, including 6 years of MODIS TERRA and 1.5 years of MODIS 
AQUA data, is available for public download through LAADS website (Level 1 and atmosphere 
archive and Distribution System, http://ladsweb.nascom.nasa.gov/data/search.html) . The 
products are accompanied by a Quality Assurance (QA) flag and color-composite RGB browse 
images for the TOA MODIS reflectance, NBRF, IBRF , cloud mask and QA. 

The algorithm has high computational efficiency. For example, full ASRVN re- 
processing of 6 years of MODIS TERRA data at 100 AERONET sites takes about 24-30 hours 
using one processor of a 2.2 GHz workstation. 

The results show very stable and reproducible NBRF and NDVI time series for any given 
pixel. The main source of errors in the developed algorithm is the variation of MODIS pixels 
with scan angle, which increase by a factor of eight from the nadir view to the edge of scan. This 
effect is partially mitigated by a 1 km resolution gridding procedure, but it is not cancelled 
entirely. In this regard, expansion of ASRVN with data from the future VIIRS instruments - 


which features constrained pixel growth — or geostationary sensors, such as the future GOES-R, 
is expected to produce a higher quality data set. Indeed, the NPOESS Land Product Validation 
Team’s initial validation plan includes an AERONET+ASRVN component. 

The ASRVN applications range from product validation and science analysis to sensor 
calibration support and long-term trending and stability studies. We believe the products from 
ASRVN fit well into the CEOS BELMANIP framework and will assist in more reliable and 
quantitative intercomparison analysis over the AERONET sites. Recently, we conducted a 
validation study of the MISR surface BRF and albedo products [10]. Because ASRVN produces 
a multi-year record for each sensor of interest, these data are useful for sensor cross-calibration 
analysis [27] and detection of long-term calibration trends. The latter are particularly important 
for climate applications. 

Our approach can be considered as an indirect validation network for current MODIS or 
MISR surface reflectance and associated products (e.g., NDVI). If supported by periodic ground 
measurements over carefully selected stable homogenous test sites with different surface 
brightness level, which would establish an absolute reference for BRF and albedo, this approach 
can also be considered a full validation that is easily expandable to global level given the 
AERONET global infrastructure. 
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APPENDIX 

A-l. Parameterized expression for the TOA radiance 

The algorithm is based on a high accuracy semi-analytical formula derived with the Green’s 
function method [23, 28], Below, r is atmospheric optical thickness, 7iS\ is spectral 
extraterrestrial solar irradiance, and s=( // = cos 6 , cp) is a vector of direction defined by zenith 
(6) and azimuthal ((p) angles. The z-axis is pointed downwards, so /Jo>0 for the solar beam and 
ju < 0 for the reflected beam. The TOA radiance L(s 0 ,s) is expressed as a sum of the atmospheric 

path radiance ( D ), and surface-reflected radiance (L s ), directly and diffusely transmitted through 
the atmosphere: 


L(s 0 ,s) = D(s 0 ,s) + L s (s 0 ,s)e /w + L d s (s 0 ,s) . 


(A-l) 


The surface-reflected radiance is written as: 


L s (s 0 ,s)^S X jUoe T/f, ° {p(s 0 ,s) + ac 0 p 1 ( j u)p 2 (ju 0 )} + — f D s (s 0 ,s')p(s',s)ju'ds' , (A-2) 

7T + 
n + 

where D s is path radiance incident on the surface, c 0 is spherical albedo of the atmosphere, and 


A (A) = J p(s',s)ds' , p 2 (p 0 ) = J P(s 0 ,s)ds . 


(A-3) 



a is a multiple reflection factor, a = (1 - q(ju 0 )c 0 ) 1 , where g is surface albedo. The diffusely 

transmitted surface-reflected radiance at the TOA is calculated from L s with the help of ID 
diffuse Green’s function of the atmosphere: 

L d s (s 0 ,s) = J G d (s 1 ,s)L s (s 0 ,s l )ds l . (A-4) 

QT 

The function nG d is often called bi-directional upward diffuse transmittance of the atmosphere. 
The method of its calculation was discussed in detail in [23]. The surface albedo is defined as a 
ratio of reflected and incident radiative fluxes at the surface: 


qW = F'*(H,)IF u " m W 


Down / 


0 / 5 


(A-5a) 


F Dow " (a,) =xS%jLI () e T,Mo + J D s (s 0 ,s')/Li'ds'= F S D "' ( a , ) + Ff' J (p 0 ) , (A-5b) 

n + 

fUp (Fo) = ^ z ju 0 e llMo q 2 (ju 0 ) + J ju'q 2 (ju')D s (s 0 ,s')ds ' , q 2 (p 0 ) = — J p(s 0 ,s)pds . (A-5c) 


These formulas give an explicit expression for the TOA radiance as a function of surface BRF. 
The accuracy of the above formulas is high, usually within a few tenths of a percent [23], Below 
we will use the TOA reflectance, which is defined as 

R a =L a /( Mo S a ). (A-6) 


A-2. Expression for the TOA reflectance using LSRT BRF model 

Based on the described semi-analytical solution, we can express TOA reflectance as an explicit 
function of parameters of the BRF model. We are using a semi-empirical Li Sparse - Ross Thick 
(LSRT) BRF model [29] as a basis. This is a linear model, represented as a sum of Lambertian, 
geometric-optical, and volume scattering components: 

p(/u 0 ,/u,<p) = k L +k G f G (iu 0 ,/u,(p) + k v f v (ii 0 ,ii,(p). (A-7) 

It uses predefined geometric functions (kernels) f G , f v to describe different angular shapes. The 
kernels are independent of the land conditions. The BRF of a pixel is characterized by a 
combination of three kernel weights, K = {k L ,k G ,k l } T . The LSRT model is used in the 
operational MODIS BRF/albedo algorithm [9]. 

The substitution of equation (A-7) into (A-l - A-5), and normalization to the reflectance units 
gives the following expressions for the surface-reflected signal (the last two terms of Eq. (A-l)): 

_ r/ 

e /f, ° {k L +k G f G ( J u 0 ,jU,(p) + k v f r (jU 0 ,jU,(p)+ac 0 p 1 (jU) p 2 (ju 0 )} 

+ ap 0 l {k L E d (a, ) + k °D g (a, ,p,<p) + k ' D' v (a ,p,(p)}, (A-8) 

R d (p 0 ,p,<p)= e^x{[k L G™{ M ) + k G Gl( Mo ,p,(p) + k v G l r {ju 0 , M ,<p)\ + 

ac 0 [, k L G av (p ) + k c G"{p) + k v G x r \p)\ p 2 {p 0 ) } 

+ «Ai ' {k ! E d (a )G av ( p ) + k G H G (a ,p,(p) + k v Hy(p 0 ,p, cp)} . (A-9) 

The surface albedo is written as: 



<j(Po) = E o (Po){Po e /Mo q 2 (jU 0 ) + k E 0 (p 0 ) + k D G (p 0 ) + k D v (ju { 0 )}. (A-10) 

Different functions of these equations represent different integrals of the incident path radiance 
( D s ) and atmospheric Green’s function (G) with the BRF kernels. They were described in [28] 
along with the numerical calculation method. Below, we give only the integral expressions: 

A ( P) = ■ k L + k' G f G ( p) + k ' fy (ju ) , 


P 2 (A ) = k L + k G f l {p f) ) + k v fy (a, ) , 
q 2 (ju 0 )=k L + k°f G (p 0 ) + k 1 fy (// 0 ), 

If.,,., 


Y 1 2ft 

p-p 0 ) =—fp'dp'f dcp'D s (p Q , //', <p' - <p 0 )f k (//', p,cp- cp r ) , 


, In 1 

Dl(pf = - J dcp’\ ju'f k (ju')D s (ju 0 ,ju'-,(p')dju ' , 

* 0 0 
0 2ft 

G av (p) = J dp x J G d {p x , p,cp - cp x )d(p x , 

-1 0 

0 2ft 

G \ X ( A> = J fl (A ) d P\ \G\ft,n,q>-q> x )d(p x , 

-1 0 

0 2ft 

G\(ii Q ,P,<P-(Po) = \dp x \G d (Hi,P,<P-<Px)f k (Po,Pi,<Px-n)d(Pi , 
-1 0 
0 2 ft 

H\(Po,P,<P-<Po) = \ d Px jG rf (A,A^-A)AUA),A,A -<Po)d<Pi ■ 


(A-l 1) 
(A- 12) 
(A- 13) 

(A- 14) 
(A- 15) 
(A- 16) 
(A- 17) 
(A- 18) 
(A- 19) 


The subscript k in the above expressions refers to either geometric-optical (G) or volumetric (V) 
kernels, and the supplementary functions of the BRF kernels are given by: 


1 i ah 

fl (n) = —\dn'\f k {n',n,q>'- <p)d<p ' , 

In 0 0 

j 0 2 k 

fl ( Po ) = — J d P\ { fk ( A > A , A - <Po ) d( p 


-1 0 
0 2 . 


^ U Aft 

fl (A') = — { P d P \f k (MfM,<P- <P')d(P ■ 

1 A 


The diffuse and total spectral surface irradiance are calculated from (A-5b) as: 

, E,(p a ) = F D -"(p l ,)l(^S i ). 


(A-20a) 

(A-20b) 

(A-20c) 

(A-21) 


Let us re-write these equations separating the kernel weights. First, single-out small terms 
proportional to the product c 0 p 2 (p 0 ) into the non-linear tenn: 


R n, (p 0 ,p)=ac 0 p 2 (p 0 ) pfp) + k L G a \ju) + £ C G“(//) + k v G x v \p ) }. 

Second, collect all remaining multiplicative factors for the kernel weights: 


(A-22) 



(A-23) 


F L (Mo,M) =(e^ + aju 0 l E d 0 (ju 0 )){e'^ + G av (p)), 

-/ - T / \ 

F k (p 0 ,/s,(p) = { e // '° f k (p Q , p,cp) + apf D l k (p Q , p,cp) } e m + 


e /m ° G l k (ju 0 ,ju,(p) + aju 0 l Hl(ju 0 ,ju,(p),k=V, G. (A-24) 

With these notations, the TOA reflectance becomes: 

R{/u 0 ,/u,(p) = R D (ju 0 ,ju,(p) + k L F L (ju 0 ,ju) + k G F G (ju 0 ,ju,(p) + k v F v \ju Q , /u,(p) + R nl (fi 0 , /u) . (A-25) 

This equation, representing TOA reflectance as an explicit function of the BRF model 
parameters, provides the means for an efficient atmospheric correction. 
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Table 1. ASRVN product suite 


Product 

Name 

Data Type 

Dimensions 

Descriptions 

MODIS 

Product 

Counterpart 

CloudMask 

DFNTUINT8 

50x50 

Cloud mask field, the 
definition of values is 
shown in Table 2. 

MOD35 cloud 
Mask 

NBRF 

DFNTFLOAT32 

50x50x7 

Bi-directional 
reflectance function 

nonnalized to SZA=45° 
and nadir view. 

MOD43B4 Nadir 
BRDF- Adjusted 
Reflectance 
(NBAR) 

Albedo 

DFNTFLOAT32 

50x50x7 

Surface albedo at a given 
solar zenith angle in 
ambient atmospheric 

conditions. 

A combination 
of MOD43B3 
black-sky and 

white-sky albedo 
weighted with 
respective 
relative direct 

and diffuse 

incident fluxes. 

IBRF 

DFNTFLOAT32 

50x50x7 

Instantaneous (or one- 
angle) BRF for specific 
viewing geometry of the 
last day of observations. 

MOD09 Surface 
reflectance 

Kiso 

DFNTFLOAT32 

50x50x7 

BRF model parameter, 
the isotropic coefficient. 

MOD43B1 

BRDF/Albedo 

Model 

Parameters 

Kvol 

DFNTFLOAT32 

50x50x7 

BRF model parameter, 
the volumetric 

coefficient. 

MOD43B1 

BRDF/Albedo 

Model 

Parameters 

Kgeo 

DFNTFLOAT32 

50x50x7 

BRF model parameter, 
the geometric-optics 

coefficient. 

MOD43B1 

BRDF/Albedo 

Model 

Parameters 

QA 

DFNT_UINT16 

50x50 

Quality assurance flags. 
The definition of each 
field can be found in 
Table 3. 

N/A 

RGB browse 
image 

N/A 

50x50 

The RGB browse image 
including TO A 

reflectance, NBRF, 

IBRF, CM and QA. 

N/A 




Table 2. Cloud mask field definition 


Value 

Definition 

1 

Clear 

2 

Possibly clear 

3 

Possibly cloudy 

4 

Cloudy 

5 

Cloud shadow 

6 

“grey” area around cloud pixel, or 
inhomogeneous aerosol. 

10 

Clear water 

11 

Clear snow 

12 

Clear Land (equivalent to 1 ) 

13 

Clear water (equivalent to 10) 




Table 3. QA field definition 


Bit field 

Bits 

Range 

Bit-code definition 

QA. overall 

0-1 

00-11 

00 - Good quality 

01 - acceptable quality 

10 - BRF parameters and NBRF are filled with 
previous results, IBRF is not created 

1 1 - No result is created 

QA. scale 

2 

0-1 

0 - Not scaled 

1 - scaled 

QA.nDelay 

3-7 

0-31 

The number of days since last update of Q-memory 

Q A. status 

8-9 

00-11 

00 - the BRF consistency parameter status =0 

01 - status =1 
10- status =2 

1 1 - status >2 (reliable solution) 

QA.adjCloud 

10 

0-1 

0 - no adjacent cloud 

1 - this pixel is adjacent to a cloudy pixel 

QA.model 

11 

0-1 

0 - the calculated IBRF is consistent with the model 
prediction 

1 - the calculated IBRF is not consistent with the 
model prediction 

Q A. cloud 

12 

0-1 

0 - the pixel is clear 

1 - the pixel is cloudy 




MODIS Tile 


Queue 



Bock/) 


Figure 1 . Structure of Queue for ASRVN processing. The Queue, designed for the sliding 
window algorithm, stores up to 16 days of gridded MODIS observations at 1 km resolution. The 
CM algorithm uses MODIS bands 1-7 and band 31, which are stored as Layers (double-indexed 
arrays) shown in the upper-left corner. A dedicated Q-memory is allocated to store the ancillary 
information for CM algorithm in Refcm structure, such as a reference clear-skies image (ref cm) 
and results of dynamic Land-Water-Snow classification ( mask_LWS ). This information is 
updated with latest measurements (day L) once given block is found cloud-free, thus adapting to 
changing surface conditions. The Q-memory also stores results of previous reliable BRF 
retrievals, or atmospheric correction (AC), for MODIS bands 1-7. 



MODISTCA 


RGB images NBRF CM IBRF AOD 



Cloud Mask: 
Reproducible 
spatial pattern 
indicates clear 
conditions 


Aerosol Filter 
detects contrails, 
thin cirrus clouds 
and non- 
homogeneous 
aerosols 



Figure 2. Example of ASRVN processing for GSFC, USA, 2006. The two left columns show the 
top-of-atmosphere (TOA) gridded RGB MODIS TERRA images for three different sequences of 
observations. The two left images are differently normalized to help distinguish clouds. The 
fourth column shows the total generated cloud mask (CM legend: clear - blue, cloud - red, 
possibly cloud - yellow, cloud shadow - dark red, aerosol fdter - grey). The last column shows 
the aerosol optical depth (AOD) retrieved from the last day of measurement, which is used in 
aerosol filter. The third and fifth columns show the results of ASRVN atmospheric correction 
(see sec. 3.3). 



For each pixel in the image: select clear-sky days 
from the queue 


Pixel level processing 





[ 

Remove outlier 


^^^nObs>3?^___^ 

Yes 

I 

BRF inversion 

Solution exists? 
JpYes 

^Max |R Cal -R Meas | 
<0.08? 

J^Yes 

Is quality 

--^^^good^ 

Yes 

X 

Increase status by 1 


Is previous 
reliable retrieval 
available? 


Calculate TOA reflectance with 
LSRT q and AERONET parameters 
for the last day (R Cal ) 


|R ca '-R Meas | <A(A)? 


Calculate scaling 
coefficient bx 


0.8<bx< 1.2? 


Is previous result more 
than 32 days old? 


Update Queue 


Compute two image-average NDVIs: 

1) NDVI q using NBRF red & NBRF Nir 

2) NDVI 1 using IBRF re d & IBRF Nir re- 
calculated for VZA=0, SZA=45° 


Is |ANDVI| >0.01 
and ARed*ANIR<0? 


Refresh Q-memory with fill value 


Image level processing 



Update queue with scaled 
LSRT parameters 



Output current 
queue results 

















Figure 3. Block-diagram of MODIS AC algorithm. 
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(a) (b) (c) 

Figure 4. . Temporal dynamics of surface NBRF for a) GSFC in year 2005; b) Mongu site in year 
2005 and c) Solar-village site in year 2006. They shown from left to right the RGB TOA MODIS 
gridded images, the RGB NBRF images and the cloud masks. The true color images are composed 
from equally weighted red, green and blue bands. . The CM legend is the same as in Figure 2. The 
red, yellow and blue dots in panel a) and c) show the location of selected “green” pixel, “urban” pixel and 
“bright” pixel, respectively. 
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Figure 4. NBRF time series for a) green pixel, GSFC site; b) urban pixel, GSFC site and c) bright 
pixel, Solar village site. R, G, B color represents red, green and blue bands, respectively; black is 
NIR band and brown is band 7 (2.1 pm). 




Figure 6. Comparison of NBRF and IBRF for GSFC, USA, day 95-122, 2006. Images shown 
from left to right are: RGB TOA MODIS gridded images, RGB NBRF images, RBG IBRF images 
and cloud mask. The CM legend is the same as in Figure 5. 
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Figure 7. NDVI time series of a) “green pixel” and b) “urban” pixel in GSFC site. 
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Figure 8. Comparison of MODIS Terra and Aqua NBRF for a bright pixel (Solar Village). Here, 
triangles and solid circles represent Aqua and Terra NBRF, respectively. Red, green and blue 
bands are shown with their respective color, black shows NIR band and brown corresponds to 
band 7 (2.1pm). 


